!
!
!===============================================================================
! module interface to FESOM2.0 for the CVMIX IDEMIX extension for the calculation 
! of hte internal wave energy and its dissipationof in Turbulent Kinetic Energy
! vertical mixing scheme --> Its based on the module interface for MPIOM programed 
! by Nils Brüggeman & Oliver Gutjahr, This module calls subroutines from the 
! CVMix library 
!
! @see Olbers D., Eden C.:
!       A Global Model for the Diapycnal Diffusivity Induced Internal Gravity Waves.
!       J. Phys. Oceanogr., 43, 1759-1779. doi: 10.1175/JPO-D-12-0207.1, 2013.
! @see Eden C., Czeschel L., Olbers D.:
!       Towards Energetically Consistent Ocean Models. 
!       J. Phys. Oceanogr., 44, 3160-3184, doi: 10.1175/JPO-D-13-0260.1, 2014.
! written by Patrick Scholz, 10.05.2019
module g_cvmix_tidal
    
    !___________________________________________________________________________
    ! module calls from cvmix library
    use cvmix_tidal,  only :  cvmix_init_tidal, cvmix_coeffs_tidal_low, cvmix_compute_Simmons_invariant_low                 
    use cvmix_kinds_and_types, only: cvmix_data_type, cvmix_global_params_type
    
    !___________________________________________________________________________
    ! module calls from FESOM
    use g_config , only: dt
    use o_param           
    use o_mesh
    use g_parsup
    use o_arrays
    use g_comm_auto 
    use g_read_other_NetCDF
    implicit none
    public
    
    !___Parameter for the init of CVMIX TIDAL___________________________________
    ! Parameters associated with the Simmons et al. scheme
    character(20) :: tidal_mixscheme="Simmons"
    
    ! units: unitless (fraction); (Gamma in Simmons et al.)
    real(kind=WP) :: tidal_efficiency=0.2                   
    
    ! units: m; zeta in Simmons et al. (used to compute the vertical deposition 
    ! function), i.e. mixing from breaking internal gravity waves, generated by 
    ! scattered barotropic tidal energy, is exponentially trapped within a 
    ! distance 'zeta' from the bottom.
    real(kind=WP) :: tidal_vert_decayscale=500.0       
    
    ! units: m^2/s; largest acceptable value for diffusivity
    real(kind=WP) :: tidal_max_coeff=50e-4            
    
    ! units: unitless (fraction);  tidal dissipation efficiency (q in Simmons 
    ! et al.), i.e. fraction of energy that dissipates locally
    real(kind=WP) :: tidal_lcl_mixfrac=0.33           
    
    ! units: m; depth of the shallowest column where tidal mixing is 
    ! computed (positive below surface)
    real(kind=WP) :: tidal_depth_cutoff=0.0              

    ! filelocation for idemix bottom forcing 
    character(200):: tidal_botforc_file = '/work/ollie/pscholz/FORCING/IDEMIX/tidal_energy_gx1v6_20090205_rgrid.nc'
    
    namelist /param_tidal/ tidal_mixscheme, tidal_efficiency, tidal_vert_decayscale, &
                           tidal_max_coeff, tidal_lcl_mixfrac,tidal_depth_cutoff,    &
                           tidal_botforc_file
    
    type(cvmix_global_params_type) :: CVmix_tidal_params ! reads the 'Prandtl' number in
    
    !___________________________________________________________________________
    ! CVMIX-TIDAL variables
    real(kind=WP), allocatable, dimension(:,:) :: tidal_Av
    real(kind=WP), allocatable, dimension(:,:) :: tidal_Kv
    real(kind=WP), allocatable, dimension(:)   :: tidal_forc_bottom_2D
    
    contains
    !
    !
    !
    !===========================================================================
    ! allocate and initialize IDEMIX variables --> call initialisation 
    ! routine from cvmix library
    subroutine init_cvmix_tidal
        
        character(len=100) :: nmlfile
        logical            :: file_exist=.False.
        integer            :: node_size
        
        !_______________________________________________________________________
        if(mype==0) then
            write(*,*) '____________________________________________________________'
            write(*,*) ' --> initialise cvmix_TIDAL'
            write(*,*)
        end if
            
        !_______________________________________________________________________
        ! allocate + initialse idemix arrays --> with size myDim_nod2D+eDim_nod2D
        node_size=myDim_nod2D+eDim_nod2D
        allocate(tidal_Av(nl,node_size))
        allocate(tidal_Kv(nl,node_size))
        tidal_Av(:,:)     = 0.0_WP
        tidal_Kv(:,:)     = 0.0_WP
        
        allocate(tidal_forc_bottom_2D(node_size))
        tidal_forc_bottom_2D(:) = 0.0_WP
        
        !_______________________________________________________________________
        ! read cvmix namelist file 
        nmlfile ='namelist.cvmix'    ! name of ocean namelist file
        ! check if cvmix namelist file exists if not use default values 
        file_exist=.False.
        inquire(file=trim(nmlfile),exist=file_exist) 
        if (file_exist) then
            open(20,file=trim(nmlfile))
                read(20,nml=param_tidal)
            close(20)
        else
            write(*,*) '     could not find namelist.cvmix, will use default values !'    
        end if    
        
        !_______________________________________________________________________
        if (mype==0) then
            write(*,*) "     tidal_mixscheme       = ", tidal_mixscheme
            write(*,*) "     tidal_efficiency      = ", tidal_efficiency
            write(*,*) "     tidal_vert_decayscale = ", tidal_vert_decayscale
            write(*,*) "     tidal_max_coeff       = ", tidal_max_coeff
            write(*,*) "     tidal_lcl_mixfrac     = ", tidal_lcl_mixfrac
            write(*,*) "     tidal_depth_cutoff    = ", tidal_depth_cutoff
            write(*,*) "     tidal_botforc_file    = ", tidal_botforc_file
            write(*,*)
        end if
        
        !_______________________________________________________________________
        ! read bottom near tidal forcing from cesm data set --> file 
        ! from N. Brüggemann interpoalted to regular grid
        file_exist=.False.
        inquire(file=trim(tidal_botforc_file),exist=file_exist) 
        if (file_exist) then
            if (mype==0) write(*,*) ' --> read TIDAL near tidal bottom forcing'
            call read_other_NetCDF(tidal_botforc_file, 'wave_dissipation', 1, tidal_forc_bottom_2D, .true.) 
            !!PS ! convert from W/m^2 to m^3/s^3
            !!PS tidal_forc_bottom_2D  = tidal_forc_bottom_2D/density_0
            ! --> the tidal energy for dissipation is divided by rho0 in 
            !     cvmix_tidal.f90 (subroutine 'cvmix_compute_Simmons_invariant_low')
            
        else
            if (mype==0) then
                write(*,*) '____________________________________________________________________'
                write(*,*) ' ERROR: TIDAL bottom forcing file not found! Cant apply TIDAL'
                write(*,*) '        vertical mixing parameterisation! '
                write(*,*) '        --> check your namelist.cvmix, tidal_botforc_file &  '
                write(*,*) '____________________________________________________________________'
            end if 
            call par_ex(0)
        end if 
        
        !_______________________________________________________________________
        ! initialise TIDAL parameters
        call cvmix_init_tidal(mix_scheme           = tidal_mixscheme,        &
                              efficiency           = tidal_efficiency,       &
                              vertical_decay_scale = tidal_vert_decayscale,  &
                              max_coefficient      = tidal_max_coeff,        &
                              local_mixing_frac    = tidal_lcl_mixfrac,      &
                              depth_cutoff         = tidal_depth_cutoff)
        
    end subroutine init_cvmix_tidal
    !
    !
    !
    !===========================================================================
    ! calculate TIDAL mixing parameterisation
    subroutine calc_cvmix_tidal
    
        integer       :: node, elem, node_size
        integer       :: nz, nln
        integer       :: elnodes(3)
        real(kind=WP) :: simmonscoeff, vertdep(nl)
        
        !_______________________________________________________________________
        node_size = myDim_nod2D
        do node = 1,node_size
            nln = nlevels_nod2D(node)-1
            vertdep = 0.0_WP
            
            !___________________________________________________________________
            ! Compute the time-invariant portion of the tidal mixing coefficient
            ! using the Simmons et al.(2004) scheme.
            call cvmix_compute_Simmons_invariant_low(         &
                 nlev            = nln ,                      &
                 energy_flux     = tidal_forc_bottom_2D(node),& !in W m-2  
                 rho             = density_0,                 &
                 SimmonsCoeff    = simmonscoeff,              &
                 VertDep         = vertdep(1:nln),                & ! vertical deposition function 
                 zw              = zbar_3d_n(1:nln+1,node),         & 
                 zt              = Z_3d_n(1:nln,node))
                
            !___________________________________________________________________
            ! Computes vertical diffusion coefficients for tidal mixing 
            ! parameterizations.
            call cvmix_coeffs_tidal_low(                      &
                 Mdiff_out       = tidal_Av(:,node),          &
                 Tdiff_out       = tidal_Kv(:,node),          &
                 Nsqr            = bvfreq(:,node),            & !FIXME: limit to N2 > 10^-8 ? as in Simmons et al.
                 OceanDepth      = -zbar_n_bot(node),         & !FIXME: neglecting free surface contribution
                 SimmonsCoeff    = simmonscoeff,              &
                 vert_dep        = vertdep(:),                &
                 nlev            = nln,                       &
                 max_nlev        = nl-1,                      &
                 CVmix_params    = CVmix_tidal_params) ! FIXME: Simmons et al. use Prandtl=10.0 (atm its 1.0)
                 
!!PS             if (mype==0) then 
!!PS                 write(*,*) 'SimmonsCoeff = ',simmonscoeff
!!PS                 write(*,*) 'tidalforcbot = ',tidal_forc_bottom_2D(node)
!!PS                 write(*,*) 'VertDep      = ',vertdep
!!PS                 write(*,*) 'bvfreq       = ',bvfreq(:,node)
!!PS                 write(*,*) 'tidal_Kv     = ',tidal_Kv(1:nln+1,node)
!!PS                 write(*,*) 'tidal_Av     = ',tidal_Av(1:nln+1,node)
!!PS             end if 
        end do !-->do node = 1,node_size
            
        !_______________________________________________________________________
        ! add tidal diffusivity to main model diffusivity Kv
        ! MPIOM note 1: simplification, other option is to explicitly calculate 
        !               the tidal scheme only below the mixed layer (as a third 
        !               option in the KPP scheme for the scheme below the mixed 
        !               layer, and then set 'depth_cutoff' to the depth level of 
        !               the bottom of the mixed layer
        ! 
        ! MPIOM note 2: background diffusivities were already added in the mixed layer 
        !               scheme (KPP)
        call exchange_nod(tidal_Kv)
        Kv = Kv + tidal_Kv
            
        !_______________________________________________________________________
        ! add tidal viscosity to main model diffusivity Av -->interpolate 
        ! therefor from nodes to elements
        call exchange_nod(tidal_Av)
        do elem=1, myDim_elem2D
            elnodes=elem2D_nodes(:,elem)
            do nz=1,nlevels(elem)-1
                Av(nz,elem) = Av(nz,elem) + sum(tidal_Av(nz,elnodes))/3.0_WP    ! (elementwise)                
            end do
        end do
        
    end subroutine calc_cvmix_tidal
end module g_cvmix_tidal
